home *** CD-ROM | disk | FTP | other *** search
/ The AGA Experience 3 / AGA Experience Volume 3 (1997)(NFA - SAdENESS)[!].iso / software / utilities / graphics / raylab / source / texture.c < prev    next >
C/C++ Source or Header  |  1996-09-01  |  19KB  |  831 lines
  1. /*
  2.     name:    texture.c
  3.  
  4.     Textures
  5.     --------
  6.  
  7.     Texture-calculation and texture-handling.
  8.  
  9.  
  10.     This source-code is part of the RayLab 1.1 package, and it is provided
  11.     for your compiling pleasure.  You may use it, change it, re-compile it
  12.     etc., as long as nobody else but you receive the changes/compilations
  13.     that you have made!
  14.  
  15.     You may not use any part(s) of this source-code for your own productions
  16.     without the permission from the author of RayLab. Please read the legal
  17.     information found in the users documentation for RayLab for more details.
  18.  
  19. */
  20.  
  21. #include  <stdlib.h>
  22.  
  23. #include  "defs.h"
  24. #include  "extern.h"
  25.  
  26.  
  27. #define ntabsize 32    /* May be too big?? */
  28.  
  29. /*  Use "global" noise table for noise functions (filled out with InitNoiseTable()) */
  30.  
  31.     unsigned short    noise_table[ntabsize][ntabsize][ntabsize];
  32.  
  33.  
  34.  
  35.  
  36. void CreateDefTexture(TEXTURE *t)
  37. {
  38.     t->CMap.Colors[0].r=1.0; t->CMap.Colors[0].g=0.3; t->CMap.Colors[0].b=0.0;
  39.     t->CMap.Colors[1].r=0.0; t->CMap.Colors[1].g=0.0; t->CMap.Colors[1].b=0.0;
  40.     t->CMap.Bounds[0]=0.0; t->CMap.Bounds[1]=1.0;
  41.     t->CMap.LastBound=1;
  42.     t->Pattern=PATTERN_NONE;
  43.     t->ImageType=IMG_NONE;
  44.     t->ImageNum=0;
  45.     t->Turbulence=0.0;
  46.     t->Reflect.r=t->Reflect.g=t->Reflect.b=0.0;
  47.     t->Filter.r=t->Filter.g=t->Filter.b=0.0;
  48.     t->Ior=1.0;
  49.     t->Ambient=0.2;
  50.     t->Diffuse=0.8;
  51.     t->Phong=0.3; t->PhongSize=10.0;
  52.     ClearTransform(&t->Transform);
  53. }
  54.  
  55.  
  56.  
  57. void CopyTexture(TEXTURE *t2, TEXTURE *t1)
  58. {
  59.     short    i,lbound;
  60.  
  61.     lbound=t1->CMap.LastBound;
  62.     for(i=0;i<=lbound;i++) {
  63.         t2->CMap.Bounds[i]=t1->CMap.Bounds[i];
  64.         CopyColor(&t2->CMap.Colors[i],&t1->CMap.Colors[i]);
  65.     }
  66.     t2->CMap.LastBound=lbound;
  67.     t2->Pattern=t1->Pattern;
  68.     t2->ImageType=t1->ImageType;
  69.     t2->ImageNum=t1->ImageNum;
  70.     t2->Turbulence=t1->Turbulence;
  71.     CopyColor(&t2->Reflect,&t1->Reflect);
  72.     CopyColor(&t2->Filter,&t1->Filter);
  73.     t2->Ior=t1->Ior;
  74.     t2->Ambient=t1->Ambient;
  75.     t2->Diffuse=t1->Diffuse;
  76.     t2->Phong=t1->Phong;
  77.     t2->PhongSize=t1->PhongSize;
  78.     CopyTransform(&t2->Transform,&t1->Transform);
  79. }
  80.  
  81.  
  82.  
  83. void CopyColor(COLOR *c2, COLOR *c1)
  84. {
  85.     c2->r=c1->r; c2->g=c1->g; c2->b=c1->b;
  86. }
  87.  
  88.  
  89. /**************************************************************
  90.  *
  91.  *     Here are some useful (more or less) texture-patterns.
  92.  *
  93.  **************************************************************/
  94.  
  95. void GetSurfaceColor(COLOR *Color, TEXTURE *t, POINT *ip)
  96. {
  97.     double    MapValue;
  98.     int    indx,foundindx;
  99.     COLORMAP *cm;
  100.     POINT    temppoint,imgpoint;
  101.  
  102.     RevTransformPoint(&temppoint,&t->Transform,ip);            /* Transform point */
  103.     if(t->Turbulence>EPSILON) {
  104.         PTurbulence(&temppoint,&temppoint,t->Turbulence);    /* Add some turbulence to the texture */
  105.     }
  106.  
  107.     Color->r=-1.0;
  108.     if(t->ImageType!=IMG_NONE) {            /* Image texture */
  109.         switch(t->ImageType) {
  110.         case IMG_24BIT:
  111.             RevTransformPoint(&imgpoint,&Img24Array[(int)t->ImageNum]->Transform,&temppoint);
  112.             GetImg24Color(Color,Img24Array[(int)t->ImageNum],&imgpoint);
  113.             break;
  114.         case IMG_8BIT:
  115.             RevTransformPoint(&imgpoint,&Img8Array[(int)t->ImageNum]->Transform,&temppoint);
  116.             GetImg8Color(Color,Img8Array[(int)t->ImageNum],&imgpoint);
  117.             break;
  118.         default:
  119.             break;
  120.         }
  121.     }
  122.     if(Color->r<0.0) {            /* Normal builtin pattern */
  123.         switch(t->Pattern) {        /* (also if outside of non-tiled image) */
  124.         case PATTERN_CHECKER:
  125.             MapValue=PatternChecker(&temppoint);
  126.             break;
  127.         case PATTERN_CIRCLES:
  128.             MapValue=PatternCircles(&temppoint);
  129.             break;
  130.         case PATTERN_RINGS:
  131.             MapValue=PatternRings(&temppoint);
  132.             break;
  133.          case PATTERN_SPOTS:
  134.             MapValue=PatternSpots(&temppoint);
  135.             break;
  136.         case PATTERN_GRADIENT:
  137.             MapValue=PatternGradient(&temppoint);
  138.             break;
  139.         case PATTERN_MARBLE:
  140.             MapValue=PatternMarble(&temppoint);
  141.             break;
  142.         case PATTERN_SOFTMARBLE:
  143.             MapValue=PatternSoftmarble(&temppoint);
  144.             break;
  145.         case PATTERN_SQUARES:
  146.             MapValue=PatternSquares(&temppoint);
  147.             break;
  148.         case PATTERN_BLURB:
  149.             MapValue=DNoise(&temppoint);
  150.             break;
  151.         case PATTERN_MANDEL:
  152.             MapValue=PatternMandel(&temppoint);
  153.             break;
  154.         case PATTERN_WOOD:
  155.             MapValue=PatternWood(&temppoint);
  156.             break;
  157.         case PATTERN_ANGULAR:
  158.             MapValue=PatternAngular(&temppoint);
  159.             break;
  160.         case PATTERN_NONE:
  161.         default:
  162.             MapValue=0.0;
  163.             break;
  164.         }
  165.  
  166.         cm=&t->CMap;
  167.         foundindx=FALSE;
  168.  
  169.         if(MapValue>EPSILON) {
  170.         indx=0;
  171.         while((indx<cm->LastBound)&&(foundindx==FALSE)) {
  172.             if((MapValue>=cm->Bounds[indx])&&(MapValue<=cm->Bounds[indx+1]))
  173.             foundindx=TRUE;
  174.             else indx++;
  175.         }
  176.         }
  177.         if(foundindx==TRUE) {
  178.         Color->r=(cm->Colors[indx].r*(cm->Bounds[indx+1]-MapValue)+cm->Colors[indx+1].r*(MapValue-cm->Bounds[indx]))/(cm->Bounds[indx+1]-cm->Bounds[indx]);
  179.         Color->g=(cm->Colors[indx].g*(cm->Bounds[indx+1]-MapValue)+cm->Colors[indx+1].g*(MapValue-cm->Bounds[indx]))/(cm->Bounds[indx+1]-cm->Bounds[indx]);
  180.         Color->b=(cm->Colors[indx].b*(cm->Bounds[indx+1]-MapValue)+cm->Colors[indx+1].b*(MapValue-cm->Bounds[indx]))/(cm->Bounds[indx+1]-cm->Bounds[indx]);
  181.         }
  182.         else {
  183.         Color->r=cm->Colors[0].r;
  184.         Color->g=cm->Colors[0].g;
  185.         Color->b=cm->Colors[0].b;
  186.         }
  187.     }
  188. }
  189.  
  190.  
  191. /**************************************************************
  192.  *
  193.  *     Builtin patterns.
  194.  *
  195.  **************************************************************/
  196.  
  197.  
  198. double PatternChecker(POINT *ip)
  199. {
  200.     return( (double)(((long)floor(ip->x+EPSILON)+(long)floor(ip->y+EPSILON)+(long)floor(ip->z+EPSILON))&1L) );
  201. }
  202.  
  203.  
  204. double PatternCircles(POINT *ip)
  205. {
  206.     return(fmod(sqrt(ip->x*ip->x+ip->y*ip->y+ip->z*ip->z),1.0));
  207. }
  208.  
  209.  
  210. double PatternRings(POINT *ip)
  211. {
  212.     return(fmod(sqrt(ip->x*ip->x+ip->y*ip->y),1.0));
  213. }
  214.  
  215.  
  216. double PatternSpots(POINT *ip)
  217. {
  218.     double    d;
  219.     POINT    dgrid;
  220.  
  221.     dgrid.x=ip->x-floor(ip->x+0.5);
  222.     dgrid.y=ip->y-floor(ip->y+0.5);
  223.     dgrid.z=ip->z-floor(ip->z+0.5);
  224.  
  225.     d=2.0*sqrt(dgrid.x*dgrid.x+dgrid.y*dgrid.y+dgrid.z*dgrid.z);
  226.     if(d>1.0) d=1.0;
  227.     return(d);
  228. }
  229.  
  230.  
  231. double PatternGradient(POINT *ip)
  232. {
  233.     return(fmod((ip->z+50000.0),1.0));
  234. }
  235.  
  236.  
  237. double PatternMarble(POINT *ip)
  238. {
  239.     double    d,dnoise,ret;
  240.     int    i;
  241.     POINT    tp;
  242.  
  243.     tp.x=ip->x*10.0;
  244.     tp.y=ip->y*5.0;
  245.     tp.z=ip->z*5.0;
  246.     dnoise=DNoise(&tp);
  247.     d=(ip->x+50000)*17.0+7.0*dnoise;
  248.     i=((int)floor(d))%17;
  249.     if(i<4) {
  250.         tp.x=ip->x*14.3;
  251.         tp.y=ip->y*20.0;
  252.         tp.z=ip->z*20.0;
  253.         ret=0.7+0.2*DNoise(&tp);
  254.     }
  255.     else if((i<9)||(i>12)) {
  256.         tp.x=ip->x*10.0;
  257.         tp.y=ip->y*10.0;
  258.         tp.z=ip->z*10.0;
  259.         d=fabs(d-floor(d*0.05882)*17-10.5)*0.1538462;
  260.         ret=0.4+0.3*d+0.2*DNoise(&tp);
  261.     } else
  262.         ret=0.2+0.2*dnoise;
  263.  
  264.     ret = (ret > 1.0 ? 1.0 : ret);        /*  0.0 <= ret <= 1.0 */
  265.  
  266.     return(ret);
  267. }
  268.  
  269.  
  270. double PatternSoftmarble(POINT *ip)
  271. {
  272.     register double    d,ret;
  273.     POINT        tp;
  274.  
  275.     tp.x = ip->x*10.0;
  276.     tp.y = ip->y*5.0;
  277.     tp.z = ip->z*5.0;
  278.     d = (ip->x+50000)*17.0+7.0*DNoise(&tp);
  279.     ret = fmod(d,17.0)*0.117647;    /* /8.5 */
  280.     if(ret>1.0) ret = 2.0-ret;    /*  0.0 <= ret <= 1.0 */
  281.  
  282.     tp.x = ip->x*(10+4*ret);
  283.     tp.y = ip->y*(5+15*ret);
  284.     tp.z = ip->z*(5+15*ret);
  285.     ret = 0.7*ret + 0.3*DNoise(&tp);
  286.  
  287.     if(ret>1.0) ret = 1.0;
  288.  
  289.     return(1.0-ret);
  290. }
  291.  
  292.  
  293. double PatternSquares(POINT *ip)
  294. {
  295.     double    d;
  296.     POINT    dgrid;
  297.  
  298.     dgrid.x=ip->x-0.5;
  299.     dgrid.y=ip->y-0.5;
  300.     dgrid.z=ip->z-0.5;
  301.     dgrid.x=2.0*fabs(dgrid.x-floor(ip->x));
  302.     dgrid.y=2.0*fabs(dgrid.y-floor(ip->y));
  303.     dgrid.z=2.0*fabs(dgrid.z-floor(ip->z));
  304.  
  305.     d=dgrid.x;
  306.     if(dgrid.y<d) d=dgrid.y;
  307.     if(dgrid.z<d) d=dgrid.z;
  308.     if(d>1.0) d=1.0;
  309.     return(1.0-d);
  310. }
  311.  
  312.  
  313. double PatternMandel(POINT *ip)
  314. {
  315.     register double    za,zb,ca,cb;
  316.     register double d,rmax;
  317.     register int    i;
  318.  
  319.     za=zb=0.0;
  320.     ca=ip->x; cb=ip->y;
  321.     rmax=4.0;
  322.  
  323.     for(i=100;i>0;i--) {
  324.         d=2.0*za*zb;
  325.         za*=za; zb*=zb;
  326.         if((za+zb)>rmax) break;
  327.         za=za-zb+ca; zb=d+cb;    /* z = z*z + c */
  328.     }
  329.  
  330.     return((double)i/(double)100.0);
  331. }
  332.  
  333.  
  334. double PatternWood(POINT *ip)
  335. {
  336.     double    r, angle;
  337.  
  338.     r=sqrt(ip->x*ip->x + ip->y*ip->y);
  339.     if(ip->y==0.0)  angle=PID2;
  340.     else angle=atan(ip->x/ip->y);
  341.     r+=0.08*sin(20.0*angle + ip->z);
  342.  
  343.     return(fmod(r,1.0));
  344. }
  345.  
  346.  
  347. double PatternAngular(POINT *ip)
  348. {
  349.     double    angle;
  350.  
  351.     if(ip->x==0.0) {
  352.         if(ip->y>0.0) angle=PID2;
  353.         else angle=-PID2;
  354.     }
  355.     else {
  356.         angle=atan(ip->y/ip->x);
  357.         if(ip->x<0.0) angle+=PI;
  358.     }
  359.     angle=fmod((angle+PIM2)*(1.0/PIM2),1.0);
  360.  
  361.     return(angle);
  362. }
  363.  
  364.  
  365. /**************************************************************
  366.  *
  367.  *     Image texture functions.
  368.  *
  369.  **************************************************************/
  370.  
  371.  
  372. void GetImg24Color(COLOR *Col,IMAGE24 *Image,POINT *ip)
  373. {
  374.     double    du,dv,x,y;
  375.     int    u,v;
  376.     register int u2,v2;
  377.     register int XSize,YSize;
  378.     COLOR    a,b,c,d;
  379.     unsigned char *Body,*Col24;
  380.  
  381.     XSize=Image->XSize; YSize=Image->YSize;
  382.     switch(Image->Maptype) {
  383.         case MAP_SPHERICAL:
  384.         MapSpherical(ip->x,ip->y,ip->z,&x,&y);
  385.         x*=(double)XSize/PIM2;
  386.         y*=(double)YSize/PI;
  387.         break;
  388.         case MAP_CYLINDRICAL:
  389.         x=MapCylindrical(ip->x,ip->y)*(double)XSize/PIM2;
  390.         if(Image->Tile==TILE_TRUE)
  391.             y=(double)YSize*(50000.0-ip->z);
  392.         else {
  393.             if((ip->z<1.0)&&(ip->z>0.0))
  394.             y=(double)YSize*(1.0-ip->z);
  395.             else
  396.             x=-1.0;
  397.         }
  398.         break;
  399.         case MAP_PLANAR:
  400.         default:
  401.         if(Image->Tile==TILE_TRUE) {
  402.             x=(double)XSize*(ip->x+50000.0); y=(double)YSize*(50000-ip->y);
  403.         }
  404.         else {
  405.             if((ip->x<1.0)&&(ip->x>0.0)&&(ip->y<1.0)&&(ip->y>0.0)) {
  406.             x=(double)XSize*(ip->x); y=(double)YSize*(1.0-ip->y);
  407.             }
  408.             else x=-1.0;
  409.         }
  410.         break;
  411.     }
  412.     if(x>=0.0) {
  413.         u=(int)floor(x); v=(int)floor(y);
  414.         du=x-(double)u; dv=y-(double)v; 
  415.         u2=(u+1)%XSize; v2=(v+1)%YSize; 
  416.         u=u%XSize; v=v%YSize;
  417.  
  418.         Body=Image->Body;
  419.         if(Image->Interpolate==INTPOL_NONE) {
  420.         Col24=&Body[3*(u+(v*XSize))];
  421.         Col->r=((double)*Col24++)*0.003921568;    /* <=> divide by 255.0 */
  422.         Col->g=((double)*Col24++)*0.003921568;
  423.         Col->b=((double)*Col24)*0.003921568;
  424.         }
  425.         else {
  426.         Col24=&Body[3*(u+(v*XSize))];
  427.         a.r=((double)*Col24++)*0.003921568;    /* <=> divide by 255.0 */
  428.         a.g=((double)*Col24++)*0.003921568;
  429.         a.b=((double)*Col24)*0.003921568;
  430.         Col24=&Body[3*(u2+(v*XSize))];
  431.         b.r=((double)*Col24++)*0.003921568;
  432.         b.g=((double)*Col24++)*0.003921568;
  433.         b.b=((double)*Col24)*0.003921568;
  434.         Col24=&Body[3*(u+(v2*XSize))];
  435.         c.r=((double)*Col24++)*0.003921568;
  436.         c.g=((double)*Col24++)*0.003921568;
  437.         c.b=((double)*Col24)*0.003921568;
  438.         Col24=&Body[3*(u2+(v2*XSize))];
  439.         d.r=((double)*Col24++)*0.003921568;
  440.         d.g=((double)*Col24++)*0.003921568;
  441.         d.b=((double)*Col24)*0.003921568;
  442.         InterpolateLinear(Col,&a,&b,&c,&d,du,dv);
  443.         }
  444.     }
  445.     else {
  446.         Col->r=-1.0;
  447.     }
  448. }
  449.  
  450.  
  451. void GetImg8Color(COLOR *Col,IMAGE8 *Image,POINT *ip)
  452. {
  453.     double    du,dv,x,y;
  454.     int    u,v;
  455.     register int u2,v2;
  456.     register int XSize,YSize;
  457.     COLOR    a,b,c,d,*CMap;
  458.     unsigned char *Body;
  459.  
  460.     XSize=Image->XSize; YSize=Image->YSize;
  461.     switch(Image->Maptype) {
  462.         case MAP_SPHERICAL:
  463.         MapSpherical(ip->x,ip->y,ip->z,&x,&y);
  464.         x*=(double)XSize/PIM2;
  465.         y*=(double)YSize/PI;
  466.         break;
  467.         case MAP_CYLINDRICAL:
  468.         x=MapCylindrical(ip->x,ip->y)*(double)XSize/PIM2;
  469.         if(Image->Tile==TILE_TRUE)
  470.             y=(double)YSize*(50000.0-ip->z);
  471.         else {
  472.             if((ip->z<1.0)&&(ip->z>0.0))
  473.             y=(double)YSize*(1.0-ip->z);
  474.             else
  475.             x=-1.0;
  476.         }
  477.         break;
  478.         case MAP_PLANAR:
  479.         default:
  480.         if(Image->Tile==TILE_TRUE) {
  481.             x=(double)XSize*(ip->x+50000.0); y=(double)YSize*(50000-ip->y);
  482.         }
  483.         else {
  484.             if((ip->x<1.0)&&(ip->x>0.0)&&(ip->y<1.0)&&(ip->y>0.0)) {
  485.             x=(double)XSize*(ip->x); y=(double)YSize*(1.0-ip->y);
  486.             }
  487.             else {
  488.             x=-1.0;
  489.             }
  490.         }
  491.         break;
  492.     }
  493.     if(x>=0.0) {
  494.         u=(int)floor(x); v=(int)floor(y);
  495.         du=x-(double)u; dv=y-(double)v; 
  496.         u2=(u+1)%XSize; v2=(v+1)%YSize; 
  497.         u=u%XSize; v=v%YSize;
  498.  
  499.         Body=Image->Body;
  500.         CMap=Image->Colormap;
  501.         if(Image->Interpolate==INTPOL_NONE) {
  502.         *Col=CMap[(int)Body[u+v*XSize]];
  503.         }
  504.         else {
  505.         a=CMap[(int)Body[u+v*XSize]];
  506.         b=CMap[(int)Body[u2+v*XSize]];
  507.         c=CMap[(int)Body[u+v2*XSize]];
  508.         d=CMap[(int)Body[u2+v2*XSize]];
  509.         InterpolateLinear(Col,&a,&b,&c,&d,du,dv);
  510.         }
  511.     }
  512.     else {
  513.         Col->r=-1.0;
  514.     }
  515. }
  516.  
  517.  
  518. void MapSpherical(double x, double y, double z, double *u, double *v)
  519. {
  520.     register double    r;
  521.  
  522.     *u=MapCylindrical(x,y);
  523.     r=sqrt(x*x+y*y+z*z);
  524.     if(r>0.0) *v=acos(z/r);
  525.     else *v=0.0;
  526. }
  527.  
  528.  
  529. double MapCylindrical(double x, double y)
  530. {
  531.     register double    r,a,b,u;
  532.     short    quadrant;
  533.  
  534.     quadrant=0;
  535.     if(x<0) {
  536.         quadrant++;
  537.         if(y<0) quadrant++;
  538.     }
  539.     else if(y<0) quadrant=3;
  540.  
  541.     r=sqrt(x*x+y*y);
  542.     if(r>0.0) {
  543.         a=fabs(x); b=fabs(y);
  544.         if(a>b) u=acos(a/r);
  545.         else    u=asin(b/r);
  546.         switch(quadrant) {
  547.         case 0:
  548.             break;
  549.         case 1:
  550.             u=PI-u;
  551.             break;
  552.         case 2:
  553.             u+=PI;
  554.             break;
  555.         case 3:
  556.             u=PIM2-u;
  557.             break;
  558.         }
  559.     }
  560.     else u=0.0;
  561.  
  562.     return(u);
  563. }
  564.  
  565.  
  566. void InterpolateLinear(COLOR *Col, COLOR *a, COLOR *b, COLOR *c, COLOR *d, double du, double dv)
  567. {
  568.     register double    dudv;
  569.  
  570.     dudv=du*dv;
  571.     Col->r = a->r + du*(b->r - a->r) + dv*(c->r - a->r) + dudv*(a->r - b->r - c->r + d->r);
  572.     Col->g = a->g + du*(b->g - a->g) + dv*(c->g - a->g) + dudv*(a->g - b->g - c->g + d->g);
  573.     Col->b = a->b + du*(b->b - a->b) + dv*(c->b - a->b) + dudv*(a->b - b->b - c->b + d->b);
  574. }
  575.  
  576.  
  577. /**************************************************************
  578.  *
  579.  *     Noise functions.
  580.  *
  581.  **************************************************************/
  582.  
  583.  
  584. void InitNoiseTable(void)
  585. {
  586.     int    x,y,z,xx,yy,zz;
  587.  
  588.     RndSeed=3426516L;    /* Ensures same table every time, on all platforms */
  589.  
  590.     for(x=0;x<ntabsize;x++) {
  591.         for(y=0;y<ntabsize;y++) {
  592.         for(z=0;z<ntabsize;z++) {
  593.             noise_table[x][y][z]=(unsigned short)floor((Jitter()+1.0)*16383.5);    /* Values range from 0 to 32767 */
  594.             if(x>=(ntabsize-1)) xx=0;
  595.             else xx=x;
  596.             if(y>=(ntabsize-1)) yy=0;
  597.             else yy=y;
  598.             if(z>=(ntabsize-1)) zz=0;
  599.             else zz=z;
  600.             noise_table[x][y][z]=noise_table[xx][yy][zz];
  601.         }
  602.         }
  603.     }
  604. }
  605.  
  606.  
  607. /*
  608.  *
  609.  *   DNoise() returns a float (double) value between 0.0 and 1.0
  610.  *
  611.  */
  612.  
  613. double DNoise(POINT *ip)
  614. {
  615.     register int  ix,iy,iz;
  616.     register unsigned short  n[8];
  617.     register double  ox,oy,oz;
  618.     double    px,py,pz;
  619.     double    n00,n01,n10,n11;
  620.     double    n0,n1,n2;
  621.  
  622.     px=ip->x-mincoord;    /* Offset x,y,z to ensure they are positive */
  623.     py=ip->y-mincoord;
  624.     pz=ip->z-mincoord;
  625.     ix=(int)px;
  626.     iy=(int)py;
  627.     iz=(int)pz;
  628.     ix=ix%(ntabsize-1);
  629.     iy=iy%(ntabsize-1);
  630.     iz=iz%(ntabsize-1);
  631.     ox=px-floor(px);
  632.     oy=py-floor(py);
  633.     oz=pz-floor(pz);
  634.  
  635.     /* interpolate to get noise value at (ix+ox,iy+oy,iz+oz)   */
  636.  
  637.     n[0]=noise_table[ix][iy][iz];
  638.     n[1]=noise_table[ix][iy][iz+1];
  639.     n[2]=noise_table[ix][iy+1][iz];
  640.     n[3]=noise_table[ix][iy+1][iz+1];
  641.     n[4]=noise_table[ix+1][iy][iz];
  642.     n[5]=noise_table[ix+1][iy][iz+1];
  643.     n[6]=noise_table[ix+1][iy+1][iz];
  644.     n[7]=noise_table[ix+1][iy+1][iz+1];
  645.  
  646.     n00=n[0]+ox*(n[4]-n[0]);
  647.     n01=n[1]+ox*(n[5]-n[1]);
  648.     n10=n[2]+ox*(n[6]-n[2]);
  649.     n11=n[3]+ox*(n[7]-n[3]);
  650.  
  651.     n0=n00+oy*(n10-n00);
  652.     n1=n01+oy*(n11-n01);
  653.     n2=(n0+oz*(n1-n0))*0.000030518;        /* 0.0 <= value <= 1.0 */
  654.  
  655. /*    n=noise_table[ix][iy][iz];
  656.     n00=n+ox*(noise_table[ix+1][iy][iz]-n);
  657.     n=noise_table[ix][iy][iz+1];
  658.     n01=n+ox*(noise_table[ix+1][iy][iz+1]-n);
  659.     n=noise_table[ix][iy+1][iz];
  660.     n10=n+ox*(noise_table[ix+1][iy+1][iz]-n);
  661.     n=noise_table[ix][iy+1][iz+1];
  662.     n11=n+ox*(noise_table[ix+1][iy+1][iz+1]-n); */
  663.  
  664.     if(n2<0.0) n2=0.0;
  665.     else if(n2>1.0) n2=1.0;
  666.  
  667.     return(n2);
  668. }
  669.     
  670.  
  671. /*
  672.  *
  673.  *   VNoise() returns a vector (v, |v| = 1) given a point p1
  674.  *
  675.  */
  676.  
  677. void VNoise(VECTOR *v, POINT *p1)
  678. {
  679.     register int    ix,iy,iz;
  680.     register unsigned short  n;
  681.     register double    ox,oy,oz;
  682.     double    px,py,pz;
  683.     double    n00,n01,n10,n11;
  684.     register double    n0,n1,f;
  685.  
  686.     f=0.000061037;
  687.  
  688.     px=p1->x+50000.0;    /* Offset x,y,z to ensure they are positive */
  689.     py=p1->y+50000.0;
  690.     pz=p1->z+50000.0;
  691.     ix=((int)floor(px))%(ntabsize-1);
  692.     iy=((int)floor(py))%(ntabsize-1);
  693.     iz=((int)floor(pz))%(ntabsize-1);
  694.     ox=px-floor(px);
  695.     oy=py-floor(py);
  696.     oz=pz-floor(pz);
  697.  
  698.     /* interpolate to get noise value at (ix+ox,iy+oy,iz+oz)   */
  699.  
  700.     n=noise_table[ix][iy][iz];
  701.     n00=n+ox*(noise_table[ix+1][iy][iz]-n);
  702.     n=noise_table[ix][iy][iz+1];
  703.     n01=n+ox*(noise_table[ix+1][iy][iz+1]-n);
  704.     n=noise_table[ix][iy+1][iz];
  705.     n10=n+ox*(noise_table[ix+1][iy+1][iz]-n);
  706.     n=noise_table[ix][iy+1][iz+1];
  707.     n11=n+ox*(noise_table[ix+1][iy+1][iz+1]-n);
  708.     n0=n00+oy*(n10-n00); n1=n01+oy*(n11-n01); v->x=1.0-(n0+oz*(n1-n0))*f;        /*  -1.0 <= v-x <= 1.0 */
  709.  
  710.     ix=(ix+5)%(ntabsize-1); iy=(iy+5)%(ntabsize-1); iz=(iz+5)%(ntabsize-1);
  711.     n=noise_table[ix][iy][iz];
  712.     n00=n+ox*(noise_table[ix+1][iy][iz]-n);
  713.     n=noise_table[ix][iy][iz+1];
  714.     n01=n+ox*(noise_table[ix+1][iy][iz+1]-n);
  715.     n=noise_table[ix][iy+1][iz];
  716.     n10=n+ox*(noise_table[ix+1][iy+1][iz]-n);
  717.     n=noise_table[ix][iy+1][iz+1];
  718.     n11=n+ox*(noise_table[ix+1][iy+1][iz+1]-n);
  719.     n0=n00+oy*(n10-n00); n1=n01+oy*(n11-n01); v->y=1.0-(n0+oz*(n1-n0))*f;
  720.  
  721.     ix=(ix+5)%(ntabsize-1); iy=(iy+5)%(ntabsize-1); iz=(iz+5)%(ntabsize-1);
  722.     n=noise_table[ix][iy][iz];
  723.     n00=n+ox*(noise_table[ix+1][iy][iz]-n);
  724.     n=noise_table[ix][iy][iz+1];
  725.     n01=n+ox*(noise_table[ix+1][iy][iz+1]-n);
  726.     n=noise_table[ix][iy+1][iz];
  727.     n10=n+ox*(noise_table[ix+1][iy+1][iz]-n);
  728.     n=noise_table[ix][iy+1][iz+1];
  729.     n11=n+ox*(noise_table[ix+1][iy+1][iz+1]-n);
  730.     n0=n00+oy*(n10-n00); n1=n01+oy*(n11-n01); v->z=1.0-(n0+oz*(n1-n0))*f;
  731.  
  732.     NormalizeVector(v,v);            /* Normalize v */
  733. }
  734.  
  735.  
  736. /*
  737.  *
  738.  *    PTurbulence() makes a random walk of 4 steps
  739.  *
  740.  */
  741.  
  742. void PTurbulence(POINT *p2, POINT *p1, double Factor)
  743. {
  744.     int    i;
  745.     double    l,o;
  746.     VECTOR    tmp1, tmp2, tmp0;
  747.  
  748.     tmp2.x=tmp2.y=tmp2.z=0.0;
  749.     VNoise(&tmp0,p1);
  750.  
  751.     l=2.0; o=0.5;
  752.  
  753.     for (i=2; i<=6; i++) {
  754.         tmp1.x=l*p1->x;
  755.         tmp1.y=l*p1->y;
  756.         tmp1.z=l*p1->z;
  757.  
  758.         VNoise(&tmp2,(POINT *)&tmp1);
  759.         tmp0.x+=o*tmp2.x;
  760.         tmp0.y+=o*tmp2.y;
  761.         tmp0.z+=o*tmp2.z;
  762.         if (i<6) {
  763.             l*=2.0;
  764.             o*=0.5;
  765.         }
  766.     }
  767.  
  768.     tmp0.x*=Factor; tmp0.y*=Factor; tmp0.z*=Factor; 
  769.     AddVector((VECTOR *)p2, (VECTOR *)p1, &tmp0);
  770. }
  771.  
  772.  
  773. /*****************************************************************
  774.  *
  775.  *    Transform point to objects original orientation
  776.  *
  777.  *****************************************************************/
  778.  
  779.  
  780. void RevTransformPoint(POINT *ReTransformedPoint, TRANSFORM *Transform, POINT *Point)
  781. {
  782.     VECTOR    tempv;
  783.     register int    i;
  784.     double    d,a,b,x,y;
  785.  
  786.     *ReTransformedPoint=*Point;
  787.     if(Transform->NumTransforms>0) {    /* Retransform the surface point */
  788.         for(i=Transform->NumTransforms-1;i>=0;i--) {
  789.         switch(Transform->Entry[i].Type) {
  790.             case TRANSFORM_SCALE:
  791.             ReTransformedPoint->x=ReTransformedPoint->x/Transform->Entry[i].Values.x;
  792.             ReTransformedPoint->y=ReTransformedPoint->y/Transform->Entry[i].Values.y;
  793.             ReTransformedPoint->z=ReTransformedPoint->z/Transform->Entry[i].Values.z;
  794.             break;
  795.             case TRANSFORM_MOVE:
  796.             ReTransformedPoint->x-=Transform->Entry[i].Values.x;
  797.             ReTransformedPoint->y-=Transform->Entry[i].Values.y;
  798.             ReTransformedPoint->z-=Transform->Entry[i].Values.z;
  799.             break;
  800.             case TRANSFORM_ROTATE:
  801.             NegVector(&tempv,&Transform->Entry[i].Values);
  802.             RevRotatePoint(ReTransformedPoint,&tempv);
  803.             break;
  804.             case TRANSFORM_WHIRL:
  805.             x=ReTransformedPoint->x;
  806.             y=ReTransformedPoint->y;
  807.             d=Transform->Entry[i].Values.x-sqrt(x*x+y*y);
  808.             if(d>0) {
  809.                 d*=-d*Transform->Entry[i].Values.y;
  810.                 a=cos(d); b=sin(d);
  811.                 ReTransformedPoint->x=a*x-b*y;
  812.                 ReTransformedPoint->y=b*x+a*y;
  813.             }
  814.             break;
  815.             case TRANSFORM_TWIST:
  816.             d=ReTransformedPoint->z*Transform->Entry[i].Values.x;
  817.             x=ReTransformedPoint->x;
  818.             y=ReTransformedPoint->y;
  819.             a=cos(d); b=sin(d);
  820.             ReTransformedPoint->x=a*x-b*y;
  821.             ReTransformedPoint->y=b*x+a*y;
  822.             break;
  823.             case TRANSFORM_NONE:
  824.             default:
  825.             break;
  826.         }
  827.         }
  828.     }
  829. }
  830.  
  831.